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Application to Large-scale System Analysis 
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Abstract —Distributed algorithms for solving coupled semidef¬ 
inite programs (SDPs) commonly require many iterations to 
converge. They also put high computational demand on the 
computational agents. In this paper we show that in case the 
coupled problem has an Inherent tree structure, it is possible 
to devise an efficient distributed algorithm for solving such 
problems. This algorithm can potentially enjoy the same effi¬ 
ciency as centralized solvers that exploit sparsity. The proposed 
algorithm relies on predictor-corrector primal-dual interior-point 
methods, where we use a message-passing algorithm to compute 
the search directions dlstrlbutedly. Message-passing here is closely 
related to dynamic programming over trees. This allows us 
to compute the exact search directions in a finite number of 
steps. Furthermore this number can be computed a priori and 
only depends on the coupling structnre of the problem. We use 
the proposed algorithm for analyzing robustness of large-scale 
uncertain systems dlstrlbutedly. We test the performance of this 
algorithm using numerical examples. 

Keywords — SDPs, distributed algorithms, primal-dual methods, 
robustness analysis, interconnected uncertain systems. 

I. Introduction 

Semidefinite programs are convex optimization problems 
that include linear matrix inequalities (LMIs) or semidefinite 
constraints. The computational complexity of solving such 
problems commonly scales badly with the number of opti¬ 
mization variables and/or the dimension of the semidehnite 
constraints in the problem. This limits our ability to solve large 
SDPs. Despite this, large SDPs are appearing more and more in 
different engineering helds, e.g., in problems related to sensor 
networks, smart grids and analysis of uncertain systems, e.g., 
see 0, 0, 0, |3T). This has been the driving force 

for devising efficient and tailored centralized solvers for such 
problems. These solvers exploit the structure in the problem 
to reduce the computational burden of solving the problem 
in a centralized manner, see e.g., 0, |Tg, gg, gg. 

Despite the success of such approaches for solving medium 
to large-scale problems, there are still problems that cannot 
be solved using centralized solvers, see e.g., 0, 0, m 
This can be due to limited available computational power 
and/or memory that prohibits us from solving the problem. 
Also it can be due to certain structural constraints, e.g., 
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privacy requirements, that obstructs us from even forming the 
centralized problem. 

For such instances, distributed algorithms may be used for 
solving the problem. These algorithms facilitate solving the 
problem using a network of computational agents, without the 
need for a centralized unit. Due to this, the computational com¬ 
plexity of these algorithms scales better, and they potentially 
enable us to address structural constraints in the problem. The 
main approach for designing distributed algorithms consists 
of two major phases. First the structure in the problem is 
exploited to decompose the problem or reformulate it as a 
coupled problem. Then first-order splitting methods are used 
for solving the resulting problem dlstrlbutedly, see e.g., 

0- This approach has been used in many applications, e.g., 
see | fT4l , gg, f^ . In the authors consider a sensor 
localization problem and use a so-called edge-based decompo¬ 
sition for reformulating the underlying SDP as a coupled one. 
They then employ alternating direction method of multipliers 
(ADMM) to solve the problem dlstrlbutedly. An optimal power 
flow problem has been considered in |jTg, where the authors 
reformulate the problem as a coupled SDP using semidehnite 
relaxation techniques. They then use ADMM to solve the 
coupled problem dlstrlbutedly. In the authors consider 
robustness analysis of large-scale interconnected uncertain 
systems. They exploit the sparsity in the interconnections 
to decompose the underlying SDP and reformulate it as a 
coupled problem. This problem is then solved dlstrlbutedly 
using algorithms that rely on proximal splitting methods. 

The algorithms designed using the aforementioned ap¬ 
proach, although effective, suffer from some issues. For in¬ 
stance, since these algorithms rely on hrst-order splitting 
methods, with convergence rates 0{l/k) or 0{l/k^) where 
k is the number of iterations, they require many iterations 
to converge to an accurate enough solution. Furthermore, 
exploiting structure and decomposing problems is commonly 
done through introduction of consensus constraints, which 
describe the coupling structure in the problem. The number 
of such constraints is commonly large for SDPs, which can 
in turn adversely affect the computational and/or convergence 
properties. Moreover the agents involved in these distributed 
algorithms need to solve an SDP at every iteration of the algo¬ 
rithm, which can potentially put a considerable computational 
burden on the agents. 

In this paper we propose a distributed algorithm for solving 
coupled SDPs with a tree structure. These SDPs are dehned 
in Section This algorithm does not suffer from any of 
the aforementioned issues. We achieve this by avoiding the 
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use of first-order splitting methods and instead rely on primal- 
dual interior-point methods, which have superior convergence 
properties. The proposed algorithm is produced by distributing 
the computations conducted at each iteration of the primal- 
dual method. Particularly, we use a message-passing algorithm 
for computing the search directions. Message passing, here, is 
closely related to non-serial dynamic programming, Q, |22| , 
p^ . We also present a similar approach for distributing the 
remaining computations at every iteration. As a consequence, 
at each iteration of the primal-dual method, the computational 
burden on each agent is very low. In fact during each iteration, 
an agent is required to factorize a relatively small matrix once 
and is required to communicate with its neighbors twelve 
times. 

The proposed algorithm in this paper is closely related to 
that of p^ . In fact, the authors in p^ use the same approach 
for devising a distributed algorithm for solving coupled non¬ 
conic problems. However, the computation of search directions 
for SDPs is not as straightforward as for non-conic problems. 
This is due to introduction of scaling matrices and their 
inverses in the KKT system, which destroys the structure in 
the problem. In order to circumvent this issue, we here put 
forth a novel way for computing the search directions at each 
iteration. This in turn enables us to use the message-passing 
algorithm for computing the search directions. 

Notice that by using this approach for computing the search 
directions, we implicitly solve the so-called augmented system. 
This is done by computing a block LDL^ factorization of 
its coefficient matrix using a fixed pivoting ordering, where 
the ordering is enforced by the coupling structure in the 
problem, p2| . This is in contrast to existing methods that 
commonly solve the so-called Schur complement system or 
normal equations. As a result, the proposed algorithm provides 
us with more stable and accurate implementation, og, m. 
Solving the augmented system is also considered in |30|, 
where the authors also compute the search directions through 
solving the augmented system by computing an LDL^ fac¬ 
torization using hxed pivoting ordering. This is particularly 
done by using regularization and iterative rehnement. In this 
paper, however, a block LDL^ factorization is computed using 
a fixed pivoting ordering without the use of regularization. 
Hence, the augmented system is solved without the need for 
iterative rehnement. 

We then use the proposed algorithm for analyzing large- 
scale interconnected uncertain systems, distributedly. This is 
made possible by exploiting the sparsity in the interconnec¬ 
tions, as outlined in Q. A similar approach was also used 
in 123) . There, the authors utilized the so-called range-space 
decomposition for reformulating the analysis problem as a 
coupled feasibility problem. They then used algorithms that 
rely on proximal splitting methods for solving it distributedly. 
We here instead use the so-called domain-space decomposition 
to reformulate the analysis problem as a coupled SDR The 
coupling structure of this coupled problem is less complicated 
than that of in j^ , and has a tree structure. This then enables 
us to use the presented distributed algorithm for solving 
the problem efficiently and distributedly. We illustrate the 
performance of the algorithm using numerical examples. 


Outline 


Next we first dehne some notations that are used throughout 
the paper. In Section |I^ we put forth a dehnition of coupled 
and loosely coupled SDPs. We review a predictor-corrector 
primal-dual interior-point method in Section III and briefly 
discuss how the structure in coupled problems is reflected in 
the computations conducted at every iteration of this method. 
Section expresses coupled problems with a tree structure 
and discusses the use of message-passing algorithm for solving 
coupled problems with a tree structure. This is then used in 
Section [V] where we present the proposed distributed algorithm 
for solving coupled SDPs with tree structure. In Section |V^ 
we discuss a decomposition approach for sparse SDPs. This 
approach is used in Section VII for reformulating the problem 
of robustness analysis of large-scale interconnected uncertain 
systems as coupled SDPs with a tree structure. We test the per¬ 
formance of the proposed distributed algorithm when applied 


to this problem using numerical experiments in Section VIII 


Finally we finish the paper with some concluding remarks in 
Section IIXI 


Notation 

We denote the set of real and complex numbers with M and 
C, and the set of mxn real and complex matrices with 
and respectively. The transpose and conjugate trans¬ 

pose of a matrix X is denoted by and X*, respectively. 
The null space of a matrix X is denoted by J^{X). With S" 
and H” we denote the set of n x n symmetric and Hermitian 
matrices. The set of integer numbers n} is denoted 

by N„. Given a set of positive integers J C N^, the matrix 
Ej G is a 0-1 matrix obtained from an n x n identity 

matrix with rows indexed by N„ \ J removed, where |J| 
denotes the number of elements in J. This means that Ejx is a 
I J|-dimensional vector that contains the elements of x indexed 
by J. We denote this vector by Xj. By and Xmn'^ we 

denote the Zth element of vector a;® and the element at row m 
and column n of matrix X® at the fcth iteration, respectively. 
Given matrices X^ for k = 1,..., N, blk diag(X^,..., 
denotes a block-diagonal matrix with blocks specified by 
the given matrices. Similarly diag(a:i,..., is a diagonal 
matrix with diagonal elements Xi,... ,X]si. Given vectors x^ 
for k = 1,... ,N, the column vector {x^, ..., x^) is all of 
the given vectors stacked. The generalized matrix inequality 
G < E[ (G < El) means that G — H is negative (semi)deflnite. 
Given a matrix X G vec(X) is an mn-dimensional 

vector that is obtained by stacking all columns of X on top 
of each other. Given two matrices X,V G X • V := 

vec(X)^ vec(V). For a symmetric matrix X G S" 

svec(X) := (Xn, V^Xai,..., y2X„i, X 22 , 

v/ 2X32, . . . , \/2Xn2, ■ ■ ■ , X„n). 

Operators mat and smat are defined as inverses of vec and 
svec, respectively. Given two matrices X and V by X 0 V 
we denote the standard Kronecker product. Given X G S", 
define U as an n(n -I-1)/2 x matrix such that U vec(X) = 
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svec(X). Then for two matrices X,Y & (gjg denotes 

the symmetrized Kronecker product that is dehned as 

x®,Y ■= (g) y + r (g) x)u'^. 

For properties of the symmetrized Kronecker product refer to 
| [M1 . Given two sets Ji and J 2 , Ji x J 2 denotes the standard 
cartesian product and by Ji x ^ Ji we denote the symmetrized 
cartesian product defined as 

J\ X s J\ .— ^ Jl ^ Jl I 3 ^ 

For these two sets Ji \ J 2 denotes the standard set minus. 
By min we denote the minimum value and with arg min we 
denote the minimizing argument of a function. By £2 we 
denote the set of n-dimensional square integrable signals, and 
TZT-C^^ represents the set of real, rational m x n transfer 
function matrices with no poles in the closed right half plane. 
A graph is denoted by Q{V,S) where V = {ui,...,u„} is 
its set of vertices or nodes and £ C V x V denotes its set of 
edges. An induced graph by y' C 1/ on Q(V,£), is a graph 
Qi(V',£') where £'= £ n (V' x V). 


II. Coupled and Loosely Coupled SDPs 


i consider a 

coupled SDP given as 

N 


minimize 

X 

i = l 

(la) 

subject to 

Q] • Xj,j_ = 6 }, j = 1,, rrii, 



i = l,...,N, 

(lb) 


A, , ^0, i = l,...,N, 

(Ic) 


where Q}, S such that [svec(( 5 j) ... svec(( 5 ^J] 
has full column rank for all i = 1 ,..., A^, with the ordered 
sets Ji C such that IJ^i Ji = j = 

with Ai G S" such that Xjk = 0 if (j, k) ^ and 

J = UiLi•= Ji Ji- This problem can be seen as 
a combination of N coupled subproblems, each of which 
dehned by the objective function kF* • X j j and constraints 
(5* • = bj for j = 1,... ,mi and Xj,j_ Y 0. Let us 

now dehne T{ij) = {k \ {i,j) G Jk}, which denotes the set 
of subproblems that are coupled in that they all depend on the 
variable Xy. Notice that agents a and b are members of 
if and only if j} C JaDJb- It is possible to provide a more 
explicit description of the coupling among the subproblems by 
decomposing Q as 

N 

minimize ^ IF% A-* ( 2 a) 

i=i 

subject to Q]»X^ = b], i = (2b) 

A“ ^ 0, i = l,...,N, (2c) 

A* = Ej.XE'^. , i = 1,..., A. (2d) 

Notice that in (|^, the objective function terms and constraints 
in (|^-(| 2 c|) are decoupled and the coupling in the problem 
is descriBed using the consensus constraints in ( |2dl i. It is 
also possible to provide a graphical representation of the 
coupling using undirected graphs. Particularly let Qs{J,£s) 


be a graph with vertex set J as dehned above and edge 
set £s = {{ii,j),iv,t)) I ^ 0}. We refer to 

this graph as the sparsity graph of the problem. Let us now 
illustrate the dehnitions above using an example given as 


minimize • A{i_ 2 , 4 }{i, 2 , 4 } + 


subject to 


Xll 

X 12 

Xl3 

0 

0 


*12 

*22 

0 

3^24 

0 


W 

A 

{1,3,4}{1.3,4} 

+ • A{4^5}{4_5)^ 

(3a) 

'*11 

*12 

*13 

0 

0 ■ 



*12 

*22 

0 

*24 

0 



*13 

0 

*33 

*34 

0 

4 0. 

(3b) 

0 

*24 

*34 

*44 

*45 



_ 0 

0 

0 

*45 

*55. 



constraint 

in (|Jbli can 

be rewritten as 


*13 

0 

0 ■ 





0 

*24 

0 





*33 

*34 

0 

= 




*34 

X44 

*45 





0 

0:45 

*55. 






E 




3:11 

2 

*12 

0 


+ E, 


*2 


*12 

*22 

*24 

'xii 

2 

*13 

0 


*24 

^44 

3 

*13 

*33 

*34 


E. 


*1 


+ E, 


*3 


*34 

3 

2M4 

3 

*45 


E. 


*2 


*45 

*55 


Ej^ ^ 0. 


with Jl = {1,2,4}, J 2 = {1,3,4} and J 3 = {4,5}. Then the 
optimal objective value of 


minimize • 

X 

A{i_2,4}{1,2,4} + 



IF^ • A{i 3 4}{i_3_4} 

+ • A{4^5}{4_5)^ 

(4a) 


xii 

2 

*12 

*14 




subject to 

X 12 

*22 

*24 

^0, 

*14 = 0, 

(4b) 


X14. 

*24 

2:44 

3 





~xii 

2 

*13 

*14 






*33 

*34 

^0, 

*14 = 0, 

(4c) 


3^14 

*34 

2:44 

3 





"0:44 

3 

*45 

^ 0, 



(4d) 


2:45 

*55 






dehnes an upperbound for the optimal objective value of 0. 
The problem in Q is a coupled SDP with smaller semidehnite 
constraints. This method for reformulating the problem is 
commonly used for cases when the original problem is either 
impossible or very difficult to solve. Notice that this problem 
is in the same format as dTJ. The sparsity graph of this problem 
is illustrated in Figure [ijwhere for instance there is an edge 
between the nodes (1,1) and (1,2) since the intersection 
between the sets X^i 1 ^ = {1, 2} and ^(1.2) = {2} is nonempty. 

In case for a coupled problem 

• \JiC\ Jj\ <^n for all i,j G N^v; 

• n I < iV for all (z, j), {v, t) G J, 

then we call this problem loosely coupled. As we will see 
later, it is possible to devise efficient distributed solvers based 
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Fig. 1. Sparsity graph for the coupled SDP in 0. 


on primal-dual interior-point methods for solving coupled and 
loosely coupled SDPs. To this end, let us first briefly review 
primal-dual interior-point methods for solving SDPs. 


111. Primal-Dual Interior-point Methods for 
Solving SDPs 

It is possible to iteratively solve a standard-form SDP, 
given as 


minimize C • X 

X 

subject to Ai • X — bi, i = 1,... ,m, (5) 

X ^ 0, 

where b G and X,Ai,C G S" such that 

[svec(^i) ... svec(Am)] has full column rank, using 
primal-dual interior-point methods. Particularly, given the it¬ 
erates >- 0, >- a primal-dual interior-point 

method generates the next iterates by 

taking a single Newton step applied to the perturbed KKT 
conditions 

Ai»X = bi, i = (6a) 

m 

ViAi + S = C, (6b) 

i = l 

XS = 51, (6c) 

together with S' 0 and X >- 0 where i5 > 0. Specifically this 
Newton step can be computed by solving the following linear 
system of equations 

A, • AX = bi - Ai • X^’‘\ i = (7a) 

m m 

AS = C-S('“)(7b) 
HoiAXS^'^'’ + X‘-^'>AS) = 51 - Hd{X^'"''S‘^^^), (7c) 

where Hd{M) = 1/2{DMD-^ + D-'^MD^), S = afi is 
the perturbation parameter with /i = X^^^ • /n denoting 
the surrogate duality gap and a G [0,1], and where is a 
modified linearization of ( |^ that ensures that the computed 
directions AS and AX are symmetric. There are different 
choices for the scaling matrix D in ( |7^ , e.g., see p8| and 
references therein. For the sake of nfevity, we limit our 


discussion to the choices presented in | [32| , | [33| , that is we 
choose D = G~^ with W = GG^ where 



This scaling is referred to as the Nesterov-Todd or NT scaling. 
In order to make the notation less complicated, from now on 
we drop the iteration index k, and we use lowercase notation 
for denoting vectorized variables or residuals, e.g., we use 
Ax as svec(AX) or r^uai as svec(i?duai)- Using symmeti'ized 
Kronecker product we can then rewrite Q more compactly as 


■ 0 

A o' 


'Av' 


Fprimal 

_1 

0 I 
U F 


Ax 

_As_ 

— 

^dual 

_ ^cent _ 


where X = [svec(Xi) ... svec{Am)]^ , U = D s D '^S, 
F = DX 08 D~^ and 


Tprimal — ^ Ax 

m 

Rau^\ = G - S ( 10 ) 
R,,,i = 5I-HD{XS), 

see p^ . One way of solving is to first solve for As as in 
As = U-i (w - C/At) , (11) 


and then solve 


1 

1 

7? 

Ax 


r 

A 0 

Av 


Fprimal 


( 12 ) 


for AX and Av, where r = rdnai — F~^r„.nt. Notice that since 
F~^U is positive definite, |38 Thm. 3.2], ( fT^ also describes 
the optimality condition for the following convex optimization 
problem 


minimize -Ax'^F ^UAx + r'^Ax 

Ax 2 (13) 

subject to AAx = rprimai- 


So it is possible to compute AX and Av by either solving 
the system of equations in ([T^ or the problem in ([10. In this 
paper we focus on predictor-corrector primal-dual methods that 
rely on modified Newton directions. In order to compute these 
directions, at each iteration, we need to solve or OD twice 
with different choices of r. We lay out a predictor-corrector 
primal-dual interior-point method in Algorithm[T] based on the 
work in fSS) . 

Remark 1: Algorithmic can detect infeasibility of the prob¬ 
lem if either the primal or dual iterates diverge. This means 
that this algorithm is unable to detect weak infeasibility, 

| |28) , and generally in such cases converges to a near-feasible 
solution, p9l . 

The major computational burden of each iteration of this 
primal-dual method concerns the computation of the predictor 
and corrector directions. Next we will investigate how the 
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Algorithm 1 Predictor-corrector Primal-dual Interior-point 
Method, 


1: Given k — 0, r € (0, 1), a € {1, 2, 3}, e > 0, efeas > 0, initial 
iterates such that ^ 0 and 0 and 

p = A<°) . S'(°V«- 

2: repeat 

Compute D. 

Predictor step: Set a — 0 and compute the search direc¬ 
tions AApred. Aupred by either solving 01 or 01 and AiS'pred 
using Ol- 

5: Compute primal and dual step sizes as 

““ (^’ A„ain((A('=))-lAAp,ed) 
ad ~ min ( 1 


A™in((SW)-iASp,ed)y ■ 


Set a = 


^ f (x(*=)-rc«pAXp„d).(S('‘)+adASp„d) A ° 

I xWtSW j ■ 

Corrector step: Having computed a compute the search 
directions AAcorr, Awcoit by either solving l|12^ or l|13[> with 


(fc) ik) / 77iik)\ — l (k) , 

^ =’'Ll-(-F^ ) ^ce„l + 

(fe)x-l 


(F 


SVGci^HO ( AApred AS'pred)), 


and AScorr using 

ASeorr = (7 


,(fe) 


svec(H-D(AAppedAVd)) - U^'°^ Axc, 


8: Compute primal and dual step sizes as above, though using 

A3^corr and AScon- 
9: Update 

+ apAX^on, 

y k+l) ^ g(k) ^ 

^(fc-l-1) ^ yh) _|_ 

10: Set fc = A: -I- 1. 

11: p = A('=) • 


12: until r 


(fe) 

primal 


■ II / (fc) \ 11 

, svec (-RdJ] 1 < efcas and fi < e. 


structure in coupled problems is reflected in ( [T3| ) and how 
this structure can be used to our advantage. 

Let us apply the primal-dual method in Algorithm H to the 
coupled SDP in Q. The perturbed KKT optimality conditions 
for this problem can be written as 


Q)»X^=b), j = (14a) 


v]Q) - smat(i;*) + = W\ 

i=i 

= 61, 

- X,,,, = 0, 


(14b) 

(14c) 

(14d) 


for z = 1,..., A^, together with 


N 


= 0, 


(15) 


and X\S^ > 0 for i Define Q* = 

[svec((5i) svec((5mi)] ■ Similar to Q, given iterates 

A and A® 0 such that they satisfy ( |2dl i, 5* 0, u* and u* 

such that Yi=iiEji ®s = 0 for alH = 1,..., A, the 

Newton step corresponding to the above system of equations 
can be computed by solving 


Q^Ay = F- Qy\ 


(16a) 


i=i 


Au'Q* — smat(Au*) -I- AS'* = 


VL* - y v]Q) + smat(z;0 - S\ (16b) 


ZA 3 

i=i 


iJi5.(AA*S* -f A*AS*) = 61 - Hd^{X^S*), 
AA* - AA. . = 0, 


(16c) 

(16d) 


for f = 1,..., A, together with 'Yi=i{Ej, ®a AF = 0, 
where the scaling matrices D* are computed as discussed 
above and in (|^, though based on the given local iterates 
A*’l^'l and S*’”. This system of equations can be rewritten 
in a more compact manner as 


0 0 

0 0 

I 

0 

0 0 


Q 0 
I -S 
0 0 
0 0 
U 0 



'Av 


^primal 


Av 


0 


Ax 

= 

I'dual 


Ax 


0 


As 


^cent 


(17) 


where Q, U and F are block-diagonal with diagonal blocks 
Q* and 


W = A* (g), 

F^ ^ 


(£i*)-^S*, 
®s (DT^, 


and = [(Aj^ (g* ... {Ej^ ®s Ej^^]. Also 

here Av, Av, Ax and As denote all the corresponding 
variables stacked, e.g., Av = (Au^,..., Au'^). Similarly 
rpiimab Tduai and Tcent denote all the primal, dual and centering 
residuals stacked, where each of the stacked terms in the 
residual vectors are based on 

yu,ai = E- Qy\ (18a) 

^duai = 1^* - > ' v]Q) + smat(A) - S\ (18b) 


Eg- 

i=i 


Rlp, = 6I-HMX^y). 


(18c) 


Similar to before, we compute the primal-dual directions by 
first solving for As as 

As* = (F*)-i (r*,„. - C/*Ax*) , 


(19) 
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or equivalently as 

As* = (F*)-i (r*,„, - W svec{ElAXEj^)) , (20) 

for i = 1,..., N. Then we solve 


■-F-iU 0 I ' 


■Ax 


r 

0 00 


Ax 


0 

Q 0 0 0 


Av 


I’primal 

I -E 0 0 


Av 


0 


where r = with r* = 

Notice that the system of equations in ( |2T] i also describes the 
optimality conditions for 

^ 1 

minimize V+ (r*)^Ax* (22a) 
Ax.Ax 2 

2=1 

subject to Q*Ax* = rp^njai, i = (22b) 

AX* - =0, i = 1,..., X. (22c) 

where (T'*)“^C/* 0 for i = 1,..., X. So the predictor and 

corrector directions can also be computed by solving ( [22l i. To 
be more precise, for the predictor directions, we solve ( |2^ , 
with tr = 0, for Axpi-ed, Axpred, Avpred and Avpred, and com¬ 
pute Aspied using ( fl^ . For the corrector directions, using the 
updated a, we compute the directions Axcon, Ax^n, Avcon- 
and Avcorr by solving ([22li with 


- {FT^ (rL. - svec(ffn‘ (AX;„,A5;„,)) (23) 

and compute Ascorr as 


As*„, = (F*)-i (r*,„- 

svec(ffn, (AX;,,A^;,J) - C/*Aai*„„), (24) 

for i = 1,...,X. As a result having computed predictor 
or corrector versions of the directions Ax, Aa;, Av and Av, 
computing ASp,.gjj and As*^,-^ can be done independently by 
X computing agents in parallel. Also notice that the coupling 
structure in ( |2^ is the same as in Q. This allows us to employ 
distributed computational algorithms to distributedly solve for 
the search directions using X collaborating agents. To illustrate 
this, note that the problem in ( |2^ can be written as 

minimize Ei(x) 

x,x ^ ^25) 

subject to Ax + Bx = c, 

with X = {Ax^,..., Ax^) and x = Ax. This problem 
can be solved distributedly using proximal splitting methods, 
e.g., ADMM, Q, pO) , p^ . The use of proximal splitting 
methods for computing the primal-dual directions has been 
considered in |^, pT) . Devising distributed algorithms for 
solving coupled SDPs that also rely on this approach can be 
seen as an extension of the use of the algorithm proposed 
in 1^ to SDPs. Even though distributed algorithms based 
on proximal splitting are effective for non-conic problems, 
they suffer from certain issues when used for solving SDPs. 
Particularly, notice that the computed search directions using 



Fig. 2. Clustered sparsity graph. 

Xl = {(1,1), (1.2), (1.4), (2,2), (2,4), (4,4)} 



X2 = {(1,1), (1,3), (1,4), (3,3), (3,4), (4,4)} K3 = {(4,4), (4,5), (5,5)) 
Fig. 3. Tree representation of the sparsity graph . 


this approach are inexact and first-order splitting methods 
generally require many iterations to compute accurate enough 
search directions. Furthermore, the number of consensus con¬ 
straints in ( |22c| i are generally large for coupled SDPs which 
can in turn adversely affect the performance and numerical 
properties of such splitting methods. Also notice that for a 
predictor-corrector primal-dual method the search directions 
are computed through solving a system of the form (22i twice. 
This means that the iterative scheme for solving (22i needs 
to be run twice at each iteration of the primal-dual method. 
Hence, distributed algorithms that rely on proximal or first- 
order splitting for computing the search directions, potentially, 
require many iterations to converge to the solution. Despite all 
such issues, in many cases such splitting methods are among 
the only resorts for distributedly solving coupled or loosely 
coupled SDPs. However for coupled problems that have an 
inherent tree structure, which is common in loosely coupled 
SDPs, we can devise an efficient algorithm for solving coupled 
SDPs. This is the focus of the upcoming sections. But first we 
express what we mean by the tree structure. 

IV. Tree Structure in Coupled Problems and 
Message Passing 

Let us reconsider the coupled SDP in (|^. Notice that for this 
problem it is possible to cluster the variables or the nodes in 
its sparsity graph as shown in Figure As can be seen from 
the figure, each of the clusters induce a complete subgraph 
on the sparsity graph. We can then provide a more compact 
representation of the sparsity graph using the tree in Figure 
Each node in this tree corresponds to each of the clusters of 
variables denoted by Ki. Furthermore, for this problem, the 
tree is such that for every two nodes i and j in the tree, KiHKj 
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is contained in all the clusters in the path connecting the two 
nodes in the tree. We refer to problems that enjoy this inherent 
structure as coupled with a tree structure. Next we lay out an 
approach for exploiting this stmcture in coupled problems. 

Let us start by describing some definitions relating to graphs. 
Consider a graph Q{V,£). A clique Ci of this graph is a 
maximal subset of V that induces a complete subgraph on 
Q, i.e., no clique is properly contained in another clique, Q. 
Assume that all cycles of length at least four of Q{V,£) have a 
chord, where a chord is an edge between two non-consecutive 
vertices in a cycle. This graph is then called chordal mi 
Ch. 4]. It is possible to make a non-chordal graph chordal by 
adding edges to the graph. The resulting graph is then referred 
to as a chordal embedding. Let Cq = {Ci,... ,Cq} denote 
the set of its cliques, where q is the number of cliques of the 
graph. Then there exists a tree defined on Cq such that for 
every Ci, Cj G Cq where i ^ j, Ci H Cj is contained in all 
the cliques in the path connecting the two cliques in the tree. 
This property is called the clique intersection property, Q, 
and trees with this property are referred to as clique trees. As 
a result it is possible to represent chordal graphs using clique 
trees. This means that in case the sparsity graph is chordal, 
it is possible to use algorithms for generating clique trees for 
chordal graphs, to extract the aforementioned tree structure 
in the problem. In fact this has been used for the coupled 
example in Q. Notice that the sparsity graph for this example 
is chordal, and the clusters marked in Figure]^ are its cliques. 
Their corresponding clique tree is depicted in Figure Also 
notice that in case the sparsity graph is not chordal, the same 
procedure can be used on its chordal embedding for extracting 
the tree structure. Coupled problems with a tree structure can 
be solved using a message-passing algorithm. Consider the 
following coupled convex optimization problem 

minimize fi(x) + f 2 (x) ~l-h fN(x), (26) 

X 

where fi : K" —>• K for j = 1,... ,N. This problem can be 
seen as a combination of N subproblems, each of which is 
defined by a term in the cost function and depends only on a 
few elements of x. Let us describe the coupling structure in 
this problem in a similar manner as we did for the coupled 
SDP in Q. That is we denote the ordered set of indices of x 
that each subproblem i depends on by J^, and we denote the 
ordered set of indices of functions that depend on xj by Ij. 
We can equivalently rewrite this problem as 

minimize fi{xjJ-\ -(27) 

where the functions fi : —>• R are lower dimensional 

descriptions of /^s such that fi{x) = fi{Ej.x) for all x G K" 
and i = 1,..., N. Let us assume that the sparsity graph of this 
problem, Qe(Ve, £s), has an inherent tree structure with a set of 
cliques Cq^ = {Ci, ..., Cg} and a clique tree, T(Vt,St)- This 
problem can be solved distributedly using the message-passing 
algorithm that utilizes the clique tree as its computational 
graph. This means that the nodes Vt = 9 } act as 

computational agents that communicate or collaborate with 
their neighbors defined by the edge set Sf In order to describe 
the message-passing algorithm, we first need to assign each 


subproblem in ( |27l l to each of the agents. We can assign a 
subproblem or function fi to an agent j if Ji C Cj. Let us 
denote the set of indices of the subproblems assigned to agent 
j by (pj. Then we can rewrite ( [26] l as 

9 

minimize E F^{x^J, (28) 


where Fi(x^J := fj The message-passing 

algorithm, much the same as dynamic programming, solves 
( |28l l by performing an upward-downward pass through the 
clique tree, see e.g., |22 Sec. 4], p6) and references therein. 
Next we show how the message-passing algorithm can be used 
for devising distributed solvers for coupled SDPs with a tree 
structure. 


V. Distributed Primal-dual Interior-point 
Methods eor Coupled SDPs 

Let us reconsider the coupled SDP in Q, and assume that 
the sparsity graph of this problem, Qs{Vs,£s), has an inherent 
tree structure with clique set Cg^ = {Ci,..., Cg} and clique 
tree T{Vt,£t). Here we propose a method that allows us to 
solve this problem distributedly over the clique tree. To this 
end, we first need to assign the constituent subproblejns of Q 
to each of the agents in the tree. Firstly define Ci C 
such that Ci Xs Ci = Cj. Then we can assign a subproblem 
i to agent j if Ji C Cj. As in Section hW let us denote 
the set of indices of subproblems assigned to agent j by 
(pj. The proposed algorithm in this section relies on primal- 
dual interior-point methods. As was discussed in Section m 
the most computationally demanding stage within the primal- 
dual method in Algorithm [T] concerns the computation of the 
predictor and corrector directions. Hence the first step for 
devising a distributed algorithm for solving coupled SDPs 
is to distribute the computation of these directions, which is 
discussed next. 


A. Distributed Computation of Primal-dual Directions Using 
Message-passing 

Recall that we can compute the predictor and corrector 
directions by solving the problem in ( [2^ for different choices 
of Firstly notice that the problem in ( |22l l is equivalent 

to the problem 


minimize -\- (C)'^Ax' (29a) 

Aa; ' ^ 2 

i=l 

subject to Q*Aa;* = rpnmai, i = 1, ■ • • ,N, (29b) 

with Ax® := svec(AA^ ^ ), that is achieved by eliminating 
the constraints in ( |22c| l. ft is then possible to compute the 
search directions by solving the problem in ( |29l ). Particularly, 
by solving this problem we compute primal variables direction 
svec(AA) and dual variables directions Av. Then we can 
construct the remaining primal and dual directions as 


AA® = AA,^,^, 

AF = (F®)-i(7®Q®Ax® -f r® - (Q®)^Au®. 


( 30 ) 




for i = I,..., N. Next theorem shows that these directions in 
fact satisfy the system of equations in (|2T|. 

Theorem 1: The primal-dual directions computed by solv¬ 
ing ( |29l ) and using ( |30l l satisfy the system of equations in (|2T|. 
Proof: Notice that any solution of ( |29l l satisfies 

(F-iU£Aa; - Q'^Av) = (31a) 

Qf AcC = Tprimal- (31b) 

By choosing AX* = AX,,j , the primal directions, Ax and 
Ax, will satisfy the third and fourth block equations in ( |2T| ). 
Furthermore, notice that by ( [3Tal l we have that 

F-^UAx - Q^Av -f r e 

So if we set 

Au* = (F*)-iC7*Q*Ax* -f r* - (Q*)^Au\ (32) 


for i = I,... ,N, not only the primal-dual iterates satisfy the 
first block equation in (|2T|, but also we have f ^Av = 0. This 
completes the proof. ■ 

Consequently, we can construct the primal-dual solutions for 
the problem in ( |22l ) by first solving the problem in ( |29l ) and 
constructing the remainder of the solution as outlined in ( |30| ). 
Notice that the coupling structure of ( [29l l is the same as that 
of Q. This means that both problems have the same sparsity 
graph and tree representation of the coupling structure. We can 
equivalently rewrite (|29|) as 


N 

minimize E (Aa;*), (33) 

where fi (Ax*) := fi (Ax*) -\-Tci (Ax*), with 

f, (Ax*) = (Ax*)'^(FO-iC/*Ax* -f (rO^Ax*, (34) 

for i = 1,..., TV, and functions Iq- for i = 1,..., N, are the 
indicator functions for the constraints in (|29b|i, i.e.. 


2'c. (Ax*) 


0 Q*Ax* = 
oo Otherwise 


This problem is in the same format as and due to its 
coupling structure, can be solved distributedly using message 
passing, see | [22| Sec. 6.2]. 

So far we have described how to distribute the computation 
of the search directions using message passing. However, it 
remains to discuss how to distributedly compute the primal and 
dual step sizes, update the perturbation parameter and decide 
on terminating the algorithm. We discuss these next. 


B. Distributed Step-size Computation and Termination Check 

The clique tree used for computing the search directions 
can also be used for performing the remaining computations 
in Algorithm [T] distributedly. Notice that the computations 
described in this section are different than that of presented 
in p2] . This is because here we rely on a predictor-corrector 
method and we are concerned with SDPs. Let us first focus 
on step size computation. Similar to the message-passing 


algorithm, in order to compute the primal and dual step sizes 
we need to perform an upward-downward pass over the clique 
tree. We start the computation from the agents at the leaves of 
the tree, where every such agent first computes 

Ap = min (^Arnin ((X^)"^AX^^^j)) , (35a) 

= min (A^in ((^^)-') , (35b) 


and communicates them to its corresponding parent. Each 
agent i that has received these quantities from their children, 
will then compute 


A], = min ( min 

iGch(i) 


in (A^) , min (Amin ((X^) ^ 

h(2) ^ 3£(i)i \ \ 


AXi 


pred 


(36a) 


Arf = mm I min 

iGch(i) 


1) (A^,),min(A„i„((5^)-iA5^„, 


(36b) 


and will communicate them to its parent. This procedure is 
then continued until we arrive at the root of the tree. At this 
point, the agent at the root computes the primal and dual step 
sizes as 

ttp := min ^1, , ad := min ^1, , (37) 

where A)) and A() are calculated as in ( [36] l. These quantities 
are then communicated downwards through the tree until they 
reach the agent at the leaves. At this point, all agents will know 
the primal and dual step sizes. So the step sizes computation 
can be done by an upward-downward pass over the tree. 
Notice that the need for computing primal and dual step sizes 
also appear in Step 7 of Algorithm [T] We can use the same 
procedure for computing the step sizes at this step by simply 
replacing the predictor directions with corrector ones. 

As can be seen from Algorithm [T] in order to compute the 
corrector directions we first need to update the parameter a 
in Step 6 of the algorithm. We can use a similar approach to 
perform this update distributedly over the clique tree. Let us 
start the computations from the leaves of the tree. Every agent 
i at the leaves will then compute and communicate 

a) = ^ (X^ + apAX^^J . {S^ + adAS^^^J, (38a) 

3&<Pi 

cr* = ^X^,S'A (38b) 

to its corresponding parent. Then every agent i that has 
received these quantities from its children computes and com¬ 
municates 

<^5 = E "*1 

iech(i) 

+ (X^' + a,AXl^J . + adASi^J, (39a) 

<^2 = ^ ^ X^ • (39b) 

iech(i) jeii>i 
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to its parent. This procedure is then continued until we reach 
the agent at the root. Then this agent also computes the 
quantities tr^ and in ( [39l l and calculates the update for 
a as 



This quantity is then communicated downwards through the 
tree until it reaches the leaves of the tree. Hence, at every 
iteration of the primal-dual method all agents will have an 
update of cr after an upward-downward pass over the clique 
tree. 

It now remains to discuss distributed computation of terms in 
the stopping criteria. This concerns the computation of primal 
and dual residuals norms together with the surrogate duality 
gap. These quantities can also be computed distributedly over 
the clique tree using an analogous approach as above. Similarly 
as before let us start the computations from the leaves of the 
tree where every such agent computes 

II 112 _ II 11 2 

rd = Ell ’"dual II ’^^ = E|| ^primal II , (41a) 

/i'= ^ X-’• S'-’ , (41b) 

jecpi 

based on the updated iterates, and communicates them to its 
parent. Then every agent i that has received the necessary 
information from its children will compute 


i 

rd = 

- E 

A 

' d 

+ E| 

1 l|2 

\r^ 

' dual ’ 

(42a) 


jech(i) 





i 

rp = 

= E 

rl 

+ E| 

1 . I|2 

\r^ 

primal ’ 

(42b) 


jech(i) 






= E 

c 

+ ^ mSf 

(42c) 


jech(i) jG<pi 


based on the updated iterates, and communicates them to the 
respective parent. This procedure is then continued until we 
reach the agent at the root, which will compute the primal and 
dual residuals as 


M 1 

' primal 

r= L 

-^+E| 

1 \\2 
\r^ 

primal ’ 

(43a) 


jech{r) 




II ^ dual 1 

i^= E 

"^.+ E| 

1 \\2 
\r^ 

^ dual ’ 

(43b) 


jech{r) 

jSpr 




C. Summary of the Algorithm and Its Computational Proper¬ 
ties 


Let us reconsider the coupled SDP in Q. Given such 
a problem and its corresponding sparsity graph, (5^(14,5^), 
we extract its tree structure based on clique set Cq^ = 
{Cl,... ,Cg}. Having done so we have the computational 
graph for our algorithm and it is possible to assign the con¬ 
stituent subproblems to each of the agents using the guidelines 
in Section IV or at the beginning of Section |V] We can now 
summarize our proposed distributed algorithm as below. 


Given fc = 0, r G (0,1), a G {1,2,3}, e > 0, Cfeas > 0, 
initial iterates 0, S’’’® 0, 

such that '^^^liEj- = Q, for 

* = 1,..., IV, and M (Ef=i 1^91) 

repeat 

for * = 1,..., g do 

Agent i, given a = 0, 

and form 

fi (Ax’) as in ( |^ for j G 4>i. 

end for 


Perform an upward-downward pass and compute the 
predictor directions using message-passing, and ( [T9] l 
and ( |^ . 

Compute the primal and dual step sizes, by performing 
an upward-downward pass through the tree as 
discussed in Section [V-BI 


Update a by performing an upward-downward pass 


through the tree as discussed in Section V-B 


for i = 1,..., g do 

Agent i reforms their subproblems with as 

in ( |2^ for j G 4>i. 

end for 


Perform an upward-downward pass and compute the 
corrector directions using message-passing, and ( [24l i 
and ( [30l l. 

Compute the primal and dual step sizes, by performing 
an upward-downward pass through the tree as 
discussed in Section IV-BI 


and the surrogate duality gap as 
1 






jgch(r) . 


(44) 


This agent will then check the stopping criteria as in Step 12 of 
Algorithm [T] If these criteria are satished, then the agent at the 
root will communicate the decision to terminate the algorithm 
downwards through the tree. Otherwise, this agent will instead 
communicate the surrogate duality gap. Agents will need this 
parameter for updating the perturbation parameter for the next 
iteration of the primal-dual method. 

So far we have expressed how to distribute the computa¬ 
tions in every iteration of the primal-dual method. Next we 
summarize the outlined distributed algorithm in this section. 


for i = 1 ,..., g do 

Agent i updates 


+ aj,AXL, 
+ adASL, 




for j 
end for 

k = k -\-l- 


Evaluate p, and the termination criteria by performing 
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an upward-downward pass through the tree and 
decide whether to terminate the algorithm. 

until the algorithm is terminated 

From the outlined algorithm, we can observe that each iteration 
of the primal-dual method is accomplished within six upward- 
downward passes through the tree. Namely, two passes for 
computing the predictor and corrector directions, two for 
computing primal and dual step sizes, one for updating cr 
and one for evaluating the stopping criteria and computing 
the surrogate duality gap. Let the height of the tree, that is 
the maximum number of edges in a path from the root to 
a leaf, be h. As a result, each iteration of the primal-dual 
method is accomplished in 6 x 2 x /i steps. Furthermore, among 
these passes the ones required for computing the predictor 
and corrector directions are by far the most computationally 
demanding ones. This is mainly because during the upward 
message-passing for these passes, every agent i needs to fac¬ 
torize a matrix, see | [2^ Sec. 6.2]. However, notice that at every 
iteration of the primal-dual method, this matrix is the same 
for the predictor and corrector directions computations. This 
means that if each agent pre-caches the factorization of this 
matrix during predictor directions computations, it can reuse 
it for corrector directions computation, see p2| Remark 8]. 
This significantly reduces the computational burden of the 
upward-downward pass for computing corrector directions. Let 
us assume that the primal-dual method converges within p 
iterations. Then the major computational burden for each agent 
concerns the computation of p factorizations of a matrix, that 
is commonly of comparatively small size for loosely coupled 
problems. This is in stark contrast to distributed algorithms 
that purely rely on first-order splitting methods, as at every 
iteration of such algorithms each agent is required to solve an 
SDR 

Remark 2: The algorithm presented in this section, can 
distributedly detect infeasibility in the sense discussed in 
Remark [T] by monitoring their local primal and dual variables. 
In case any agent detects divergence of these variables, it can 
then communicate the occurrence through the tree to terminate 
the algorithm. 

Next we discuss a class of sparse SDRs, that appear in robust¬ 
ness analysis of large-scale interconnected uncertain systems, 
and we will describe how such problems can be reformulated 
as coupled SDRs with an inherent tree structure. 

VI. Chordal Sparsity and Domain-space 
Decomposition 

In order to describe sparsity in SDRs, we first briefly discuss 
the use of graphs for expressing sparsity patterns of symmetric 
matrices. 

A. Sparsity and Semidefinite Matrices 

Consider a symmetric matrix AT £ S", and an undirected 
graph H(y,£) with V = n} and £ = {{i,j) G 

{V X V) I Xij 7 ^ 0, i j}. We refer to this graph as 
the sparsity pattern graph of X. It is also possible to use 
undirected graphs to describe partial symmetric matrices. A 


partial symmetric matrix is a symmetric matrix where only a 
subset of its elements are specified and the rest are free. For 
the symmetric matrix X this structure can be expressed using 
H{V, £) with V = {1,. .., n} and £ C [V x V). Rarticularly, 
the edge set is such that we can express the set of indices of 
specified elements using Ig = f U {(t, i) | i = 1,, n}. We 
denote the set of partial symmetric matrices over H{V,£) by 
S^. A matrix X £ is then said to be positive semidefinite 
completable if by choosing its free elements, i.e., elements 
with indices in 1/ = (V x V) \ Ig, it is possible to produce 
a positive semidefinite matrix. Such matrices play a central 
role in the upcoming discussions. Let us review a fundamental 
result concernin g se midefinite completable matrices. 

Theorem 2: |^ , Thm. 7] Let H{V,£) be a chordal graph 
with clique set = {Ci,... ,Ci\ such that clique inter¬ 
section property holds. Then X £ is positive semidefinite 
completable, if and only if 

^ 0 , i = l,...,l. ( 45 ) 

We will next discuss how this theorem can be used for 
reformulating sparse SDRs. 


B. Domain-space Decomposition 

Consider the following inequality-form SDR 


minimize c y 
y 


(46a) 


subject to E ElQ^Ej^y, + E E'^M^Ej^ a 0 (46b) 


2 = 1 


2 = 1 


where y £ K®, Q^,M^ £ and Ji C N„ for i = 

1,... ,g. Let us denote the sparsity pattern graph for the matrix 
Si=i with H{V, £). Assume that this graph is chordal, 

or that we can produce a chordal embedding by adding a few 
edges, with clique set Ch = {Ci,... ,Ci}. The dual problem 
for (|46ll is given as 


g 


minimize 

z 

2 = 1 


(47a) 

subject to 

• Q* = -Ci, i = 1,.. 

■,9, 

(47b) 


ZAO. 


(47c) 


We can observe that the only elements that affect the equality 
constraints and the cost function are the ones specified by Ig. 
The rest are only used in the semidefinite constraint. This in 


turn implies that Z 

£ S^, and 

using Theorem]^ 

allows us to 

equivalently rewrite 

(|47|) as 


minimize 

-E^a. 


(48a) 

subject to 

^AA • Q" 

= C2, i = 1,.. 

.,g, (48b) 


Ze e- ^ 

i = 1,..., 

(48c) 


This method of reformulating (|47] i as (|48] l is referred to as the 
domain-space decomposition |T^, m Notice that for every 
Ji there exists a Cj such that Ji C Cj. This is because every 



11 


set Ji induces a complete subgraphs on H{V,£), and hence 
based on the definitions of cliques, it is either a subset of a 
clique or a clique itself. Le^ us denote the set of indices of 
sets Ji that are a subset of Cj by (j)j. We can then group the 


equality constraints in (|48b| and rewrite the problem in 

(j^ as 

minimize 


(49a) 



subject to 

Zj.j- •Q^ = -Cj, j Gfi, i = 1, 

...,l, 

(49b) 


c ^ i = l,...,l. 

(49c) 


which is in the same format as 0. This problem comprises 
I subproblems. Furthermore, due to its construction has a 
chordal sparsity graph with q cliques and a clique tree that has 
the same structure as the clique ttee for iF {V,£), where instead 
of Ci, the cliques are given as Ci XgCi. In fact the chordality 
of the sparsity graph follows, since the ordering defined by 
the clique tree is also a perfect elimination ordering for this 
graph, see GZl for more details. 

Remark 3: Notice that the discussion in this section also 
extends to matrices in positive semidefinite Hermitian cones, 
m- This means that the decomposition scheme described 
here, can also be used for problems with complex data matri¬ 
ces. 

Next we will discuss robustness analysis of interconnected 
uncertain systems and will show how the approach described 
here can be used for reformulating this problem as a coupled 
SDR 


where v and A{v) are the Fourier transforms of the signals 
The uncertain system is then said to be robustly 
stable if the interconnection between G and A remains stable 
for all A G IQC(n). This can be established using the 
following theorem. 

Theorem 3: The uncertain system in ( |50l l is robustly sta¬ 
ble, if 

1) for all T G [0,1] the interconnection described in ( |50l l, 
with tA, is well-posed; 

2) for all T G [0,1], tA e IQC(n); 

3) there exists e > 0 such that 


/ 


nO'w) 


Gijoj) 

I 


:< —el, Vw G [0, ooj. 


(53) 


Proof: See f^]. ■ 

Satisfaction of the conditions in this theorem is a sufficient 
condition for robustness of the uncertain system. As a result, 
for robustness analysis of this system it is required to find 
a multiplier 11 such that A G IQC(n) and that it satisfies 
the semi-infinite LMI in ( |53l l. The condition A G IQC(n) 
commonly imposes structural constraints on 11, and hence the 
analysis problem is then to find 11 with a particular structure 
such that it satisfies ( [5^. I t is possible to do this using 
either the KYP lemma, |20|, | |35| , or approximately using 
frequency-gridding, which establishes satisfaction of over 
a finite frequencies. We utilize the latter approach later as 
it preserves the structure in the problem. Next we describe 
how this framework can be used for analyzing interconnected 
uncertain systems. 


VII. Robustness Analysis of Interconnected 
Uncertain Systems 

In this section, we discuss robustness analysis of intercon¬ 
nected uncertain systems using integral quadratic constraints 
(IQCs). We start this discussion by first reviewing the IQC 
analysis framework. 


A. Robustness Analysis using IQCs 

Consider the following uncertain system 


p = Gq, q = A(p), (50) 


where G G is the system transfer function matrix, 

and A : K"* —is a bounded and causal operator 
representing the uncertainty in the system. We can characterize 
the uncertainty in the system using IQCs. Particularly it is said 
that A satisfies the IQC defined by 11, i.e., A G IQC(n), if 



V 

A{v) 


n 


V 

A{v) 


dt > 0, 


WvGCi, 


(51) 


where If is a bounded and self-adjoint operator. This constraint 
can also be written in frequency domain as 


POO 

v{ju}) 

* 

v{juj) 

J —oo 


n(jw) 



doj > 0, 


(52) 


B. Robustness Analysis of Interconnected Uncertain Systems 
using IQCs 

An interconnected uncertain system can be viewed as a 
network of N uncertain subsystems. We describe each of these 
subsystems as 


f = Gi^f + g;^w^ 

= Gl^q^ + (54) 

f = A\p% 


where Gj,, G G^,^ G , G*, G 

G*„ G and A* : ^ It is possible to 

describe the interconnection among the subsystems using a 
0-1 matrix F as 


wA 

w'^ 


Til ri2 

r2i r22 

Cin' 

C 2 N 


rzh 



.bivi riv2 ■ 

• ^NN. 


.Z^. 


r 


(55) 


where each Fij describes which components of is connected 
to which components of w*. Let us define p = (p^,... ,p^), 
q = [q^,..., q^), w = (w^, ..., w^) and z = (z^,..z^). 
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Then we can compactly describe the entire interconnected 
uncertain system as 


w 


P — GpqQ Gpw 
^ = GzqQ + Gzw'^ 
q = A{p) 
w = Tz, 


(56) 


where G*. = diag(Gl,,..., G^) and A = 

diag(A^,..., A^). Let us assume that the interconnected 
system is nominally or internally stable, i.e., 
(/ — LG^uj)”^ G TZT-C^'^ with rh — was 

then shown in Q that the system is robustly stable if there 
exist 


n = 


nil 

n 2 i 


ni 2 

n 22 


with 

IQC 


n*. 


Tl] 


TT® TT® 

iiol iio 


such that 


diag(ni,,...,nf,) and A* 
and a diagonal matrix X >- 


G 

0 


Gpq Gpw 

* 

Rii ni2 


Gpq Gpw 

I 0 


1121 1122 


I 0 




X[-TG,, I-rG,^]<-eI. (57) 


It is possible to rewrite this problem in the following standard 
form 



Fig. 4. A chain of N uncertain subsystem. 




Fig. 5. Convergence behavior of the algorithm for analysis of a chain 
of uncertain systems. The figure on the left shows the sum of primal 
and dual residuals and the figure on the right depicts the surrogate 
duality gap. 


find y (58a) 

771 

subject to j/iQ* + TL ^ 0 (58b) 

i=l 

where Q* G for alH = 1,..., m and W G with 

d — We can equivalently rewrite the problem in 

as below 


find 


subject to 


y 

771 

1=1 


Re(Q*) -Im(Q*) 
lm(0*) Re(Q0 


+ 


W 

0 


0 

w 


(59a) 


^ 0 


(59b) 


where all the data matrices are real, 01). In case L is sparse, 
then this SDP is also sparse and can be written in the same 
format as in ([4^ with c = 0. As a result we can use 


the approach presented in Section VI for reformulating this 
problem as a coupled problem, and employ the algorithm 
presented in Section |V] for solving it. 

Remark 4: As was discussed in Remark ^the decomposi¬ 
tion can be conducted directly on 0 or (|58] l. However, we 
here choose to reformulate the problem in ( |59| l, with real data 
matrices, for ease of notation and ease of use of the algorithm 
described in Section 0 

Next we illustrate this approach and study the performance 
of the algorithm using numerical experiments. 


VIII. Numerical Experiments 

In this section we consider two examples, namely a chain 
of uncertain systems and an interconnected uncertain system 
over a so-called scale-free network. These examples are taken 
from Q. Let us start with the analysis of a chain of uncertain 
systems, as illustrated in Figure As can be seen from the 
figure, for subsystems 1 < i < N, z\w^ G and for 
subsystems i = 1,N, z'^,w'^ G K. The uncertainty in each 
subsystem i is represented using (5®, which is assumed to 
be an unknown gain in the normalized interval [—1,1]. We 
can hence describe the uncertainties as (5* G IQC(n*) with 

W = , and r,(ycc) > 0, |31|. The inter¬ 

connection matrix for this interconnected system is described 
by the nonzero blocks T^ i-i = ^ for 

where ^ ^ 


- 


f = 2,...,iV, 
= 3,...,A - 1, 


— — (Oil)- 


0 0 

and r 2 i = r ^2 = (IjO)) r^v-i.AT 
We considered the analysis problem for this system with 
N = 100 subsystems in the chain, at a single frequency 
a; = 1 radjs. We solved 10 instances of this problem with 
different transfer function matrices for the subsystems. The 
transfer function matrices for each instance were randomly 
generated using the approach presented in Q. This guar¬ 
antees that the interconnected system is robustly stable for 
all instances. Furthermore, for this problem the multiplier X 
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was chosen as X = diag(xi,... ,X 2 N- 2 )- This resulted in a 
problem in the same format as in ( |58] l, with m = 298 and 
W G S298. 

Forming ( |59l l for this analysis problem, resulted in an LMI 
with a chordal sparsity pattern, with 198 cliques where the 
largest clique was of size 8 . The clique tree over these cliques 
had a height of 99. In order to establish chordality of the 
sparsity pattern graph and generate its cliques a greedy search 
algorithm with min degree criterion was used, fl?) . If we 
now form the problem in ( |49] l, this problem will comprise 
198 subproblems and can be solved distributedly over the 
clique tree. The parameters within the primal-dual method 
were chosen to be the same for all instances and are chosen as 
a = 1, T = 0.98, e = Cfeas = 10-^2^ = 0 and = 0 

for all i = 1,..., TV, and and for i = 1,..., TV 

were chosen to be diagonal matrices with positive diagonal 
entries generated randomly with a uniform distribution in the 
interval (0.1,2). In the worst case the primal-dual method 
converged after 12 iterations. The convergence behavior of this 
instance is illustrated in Figure and as can be seen mimics 
that of a standard primal-dual method, i.e., convergence within 
10 to 50 iterations with a quadratic convergence phase, El- 
Considering the height of the tree, this algorithm then, in the 
worst case, converged after 6 x 2 x 99 x 12 = 14256 steps. 
During the run of the algorithm, each agent was required to 
compute a factorization 12 times and needed to communicate 
with its neighbors 144 times. The computations in the remain¬ 
ing steps were trivial. 

We further tested the performance of the algorithm on 
a larger example with a more complicated interconnection 
description. Particularly we used the same scale-free net¬ 
work as in Sec. 5.2] for describing the interconnections 
among the subsystems. This resulted in an extremely sparse 
interconnection matrix. The transfer function matrices for the 
subsystems were also generated using the approach presented 
in 0 . Forming ( |59| ) for this analysis problem resulted in an 
LMI that is sparse with m = 1498 and W G The 

chordal embedding for the sparsity pattern graph of this LMI 
was generated by introducing 2.4% hll-in, also using a greedy 
search algorithm, with 579 cliques. The largest of these cliques 
had a size of 261. The corresponding clique tree for this 
problem was of height 35. This means that the corresponding 
problem in ( |49l ) will comprise of 579 subproblems and can 
be solved distributedly over this clique tree. We tested the 
performance of the proposed algorithm over 10 instances of 
this problem. The parameters of the primal-dual method were 
chosen to be the same as above. In the worst case the algorithm 
converged after 14 iterations. The convergence behavior of this 
instance is illustrated in Figure As a result, in the worst 
case, the algorithm converged after 6 x 2 x 35 x 14 = 5880 
steps. During the run of the algorithm, each agent needed to 
compute a factorization only 14 times and were required to 
communicate with its neighbors 168 times. 

IX. Conclusions 

In this paper we put forth a distributed algorithm for solving 
coupled SDPs with a tree structure. The proposed algorithm, 
unlike the existing ones, does not use hrst-order splitting 




Fig. 6. Convergence behavior of the algorithm for analysis of an 
interconnected system over a scale-free network. The figure on the 
left shows the sum of primal and dual residuals and the figure on the 
right depicts the surrogate duality gap. 


methods but instead uses primal-dual interior-point methods. 
Particularly, this algorithm utilizes the inherent tree structure 
in the problem as its computational graph, and distributes 
the computations at each iteration of the primal-dual method 
among the computational agents. In order to compute the 
search directions at every iteration, we employ a message¬ 
passing algorithm. This enables us to compute the exact search 
directions in a hnite number of iterations. Furthermore, we 
showed that this number can be computed a priori and only 
depends on the height of the tree. We applied the proposed 
algorithm for solving robustness analysis of large-scale inter¬ 
connected uncertain systems, and illustrated the performance 
of the algorithm using numerical experiments. 

As was discussed in the introduction, designing distributed 
algorithms are commonly conducted in two phases. Namely, 
a decomposition or reformulation phase and a splitting phase. 
In this paper, we mainly focused on the second phase of this 
procedure, that is design of efficient methods to distribute 
the computations of solving a given coupled SDR However, 
it is possible to further improve the computational and/or 
implementation properties of the devised algorithm, by using 
the available flexibilities in decomposition or reformulation 
phase. We will explore such possibilities as future line of 
research. This will mainly concern devising heuristics for 
clique or cluster merging to reduce the overall computational 
cost of the algorithm and/or to better represent the intuitive 
properties of the problem, such as physical structure in the 
problem. 
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